    ! mix metropolis algorithm for ising+XY model on the square lattice
    ! H=\sum -Jcos(\theta_i-\theta_i+x)-J\sigma_i\sigma_i+ycos(\theta_i-\theta_i+y)   
    !J=Ja;Jb=0	
    !--------------------------------------------------------------------------
    module global
    !模型定义
    integer :: n,Nsite
    double precision :: J,pi
    integer,allocatable :: site_ising(:,:)
    integer,allocatable :: nbor(:,:)
    double precision,allocatable :: site_xy(:,:)

    !模型参数
    integer :: T_number
    double precision :: T,beta
    double precision :: T_final,T_initial,T_interval
    integer :: Thermal,sampling,interval

    !模型数据
    double precision :: energy,energy_mean,arm,asm,susp,ase,ave,speh,tbr,M_2,M_4,rem,rem_2,sqm,var
    double precision,allocatable :: energy_sampling(:),var_sampling(:)


    !随机数
    integer :: is, ns, clock
    integer, dimension(:), allocatable :: seed

    contains
    !input---------------------------------
    subroutine input
    allocate(site_ising(n,n))
    allocate(site_xy(n,n))
    allocate(nbor(Nsite,4))
    allocate(energy_sampling(sampling))
	allocate(var_sampling(sampling))
    end subroutine input
    !generate ising model------------------
    subroutine generate_ising
    do i=1,n
        do j=1,n
            site_ising(i,j)=-1;
        enddo
    enddo
    end subroutine generate_ising
    !generate xy model---------------------
    subroutine generate_xy
    do i=1,n
        do j=1,n
            call random_number(rn)
            site_xy(i,j)=0;
        enddo
    enddo
    end subroutine generate_xy
    !generate table of near neighbor-------
    subroutine neighbor
    integer :: iy,ix,ixp,ixm,iyp,iym
    do ispin=1,Nsite
        iy=((ispin-1)/n)+1;
        ix=ispin-(iy-1)*n;
        ixp=ix+1-(ix/n)*n;
        iyp=iy+1-(iy/n)*n;
        ixm=ix-1+((n-ix+1)/n)*n;
        iym=iy-1+((n-iy+1)/n)*n;
        nbor(ispin,1)=(iy-1)*n+ixp;!右邻居
        nbor(ispin,4)=(iyp-1)*n+ix;!下邻居
        nbor(ispin,3)=(iy-1)*n+ixm;!左邻居
        nbor(ispin,2)=(iym-1)*n+ix;!上邻居
    enddo
    end subroutine neighbor
	
    !flip spin-----------------------------
    subroutine flip_spin
    double precision :: E_new,E_old,deltaE,alpha,temp1
    double precision :: t1,t2,t3,t4,t5,t6,t7,t8
    integer :: current_site,right,up,left,down
    !pick a site
    call random_number(rn)
    current_site=1+floor((Nsite)*rn); iy=((current_site-1)/n)+1;   ix=current_site-(iy-1)*n;
    right=nbor(current_site,1);       iy1=((right-1)/n)+1;         ix1=right-(iy1-1)*n;
    up=nbor(current_site,2);          iy2=((up-1)/n)+1;            ix2=up-(iy2-1)*n;
    left=nbor(current_site,3);        iy3=((left-1)/n)+1;          ix3=left-(iy3-1)*n;
    down=nbor(current_site,4);        iy4=((down-1)/n)+1;          ix4=down-(iy4-1)*n;
    call random_number(rn)
    alpha=rn*2*pi;
    ! H=\sum -J(1+\sigma_i\sigma_j)cos(\theta_i-\theta_j)
	! H=\sum -Jcos(\theta_i-\theta_i+x)-J\sigma_i\sigma_i+ycos(\theta_i-\theta_i+y) 
    t1=-1*J*cos(site_xy(iy,ix)-site_xy(iy1,ix1));
    
    t2=-J*site_ising(iy,ix)*site_ising(iy2,ix2))*cos(site_xy(iy,ix)-site_xy(iy2,ix2));

    t3=-1*J*cos(site_xy(iy,ix)-site_xy(iy3,ix3));
    
    t4=-J*site_ising(iy,ix)*site_ising(iy4,ix4))*cos(site_xy(iy,ix)-site_xy(iy4,ix4));
    
	E_old=t1+t2+t3+t4 ;

    temp1=2*alpha-site_xy(iy,ix); 
    if (temp1>2*pi)then
        temp1=temp1-2*pi;
    endif
    if (temp1<0)then
        temp1=temp1+2*pi;
    endif
    t1=-1*J*cos(temp1-site_xy(iy1,ix1));
    
    t2=-J*site_ising(iy,ix)*site_ising(iy2,ix2))*cos(temp1-site_xy(iy2,ix2));

    t3=-1*J*cos(temp1-site_xy(iy3,ix3));
    
    t4=-J*site_ising(iy,ix)*site_ising(iy4,ix4))*cos(temp1-site_xy(iy4,ix4));
    
	E_new=t1+t2+t3+t4 ;

    deltaE=E_new-E_old;
    call random_number(rn)
    if (rn<exp(-1*deltaE*beta))then
        site_xy(iy,ix)=temp1;
    endif
    !--------------------------------------------------------------------------
    !pick a site
    call random_number(rn)
    current_site=1+floor((Nsite)*rn); iy=((current_site-1)/n)+1;   ix=current_site-(iy-1)*n;
    right=nbor(current_site,1);       iy1=((right-1)/n)+1;         ix1=right-(iy1-1)*n;
    up=nbor(current_site,2);          iy2=((up-1)/n)+1;            ix2=up-(iy2-1)*n;
    left=nbor(current_site,3);        iy3=((left-1)/n)+1;          ix3=left-(iy3-1)*n;
    down=nbor(current_site,4);        iy4=((down-1)/n)+1;          ix4=down-(iy4-1)*n;
    ! H=\sum -J(1+\sigma_i\sigma_j)cos(\theta_i-\theta_j)
    t1=-1*J*cos(site_xy(iy,ix)-site_xy(iy1,ix1));
    
    t2=-J*site_ising(iy,ix)*site_ising(iy2,ix2))*cos(site_xy(iy,ix)-site_xy(iy2,ix2));

    t3=-1*J*cos(site_xy(iy,ix)-site_xy(iy3,ix3));
    
    t4=-J*site_ising(iy,ix)*site_ising(iy4,ix4))*cos(site_xy(iy,ix)-site_xy(iy4,ix4));
    
	E_old=t1+t2+t3+t4 ;

  

    t1=-1*J*cos(site_xy(iy,ix)-site_xy(iy1,ix1));
    
    t2=J*site_ising(iy,ix)*site_ising(iy2,ix2))*cos(site_xy(iy,ix)-site_xy(iy2,ix2));

    t3=-1*J*cos(site_xy(iy,ix)-site_xy(iy3,ix3));
    
    t4=J*site_ising(iy,ix)*site_ising(iy4,ix4))*cos(site_xy(iy,ix)-site_xy(iy4,ix4));
    
	E_new=t1+t2+t3+t4 ;

   
    deltaE=E_new-E_old;
    call random_number(rn)
    if (rn<exp(-1*deltaE*beta))then
        site_ising(iy,ix)=-site_ising(iy,ix);
    endif
    end subroutine flip_spin
    !calculate-----------------------------
    subroutine calculate
    double precision :: s_x1,s_y1,s_x2,s_y2
    !H=\sum -J(1+\sigma_i\sigma_j)cos(\theta_i-\theta_j)
    energy=0.0;
	m1=0.0	!updete9.24
	!rho_s=0.0   !update9.26
	var=0.0;
    do j=1,n
        do i=1,n
            do k=1,2
                next_site=nbor((j-1)*n+i,k);   iy=((next_site-1)/n)+1;   ix=next_site-(iy-1)*n;
                t1=-1*J*cos(site_xy(iy,ix)-site_xy(iy1,ix1));
    
                t2=-J*site_ising(iy,ix)*site_ising(iy2,ix2))*cos(site_xy(iy,ix)-site_xy(iy2,ix2));

                t3=-1*J*cos(site_xy(iy,ix)-site_xy(iy3,ix3));
    
                t4=-J*site_ising(iy,ix)*site_ising(iy4,ix4))*cos(site_xy(iy,ix)-site_xy(iy4,ix4));
    
	            energy=t1+t2+t3+t4 ;
                m1=m1+site_ising(j,i)		!update9.24		
            enddo
			!计算超流密度
			current_site=(j-1)*n+i;
			right=nbor(current_site,1);       iy1=((right-1)/n)+1;         ix1=right-(iy1-1)*n;
            up=nbor(current_site,2);          iy2=((up-1)/n)+1;            ix2=up-(iy2-1)*n;
            left=nbor(current_site,3);        iy3=((left-1)/n)+1;          ix3=left-(iy3-1)*n;
            down=nbor(current_site,4);        iy4=((down-1)/n)+1;          ix4=down-(iy4-1)*n;
		    var=var+sin(site_xy(j,i)-site_xy(iy3,ix3))!+cos(pi/4)*sin(site_xy(j,i)-site_xy(iy2,ix2))
        enddo
    enddo
	energy=energy/Nsite;
	rem=abs(float(m1)/Nsite)! absolute magnetization per site,updete9.24
	rem_2=float(m1)/Nsite   !update9.15求m,updete9.24
    var=var*var/Nsite;
	
	
	
	
	
    end subroutine calculate
	
	
    !generate random number----------------
    subroutine init_random_seed()
    call random_seed(size = ns)
    allocate(seed(ns))
    call system_clock(count=clock)
    seed = clock + 37 * (/ (is - 1, is = 1, ns) /)
    call random_seed(put = seed)
    deallocate(seed)
    end subroutine init_random_seed
    
    end module global

    program mcmc_mix_metropolis
    use global
    implicit none
    integer :: i,ii,iii

    !读取:维度格点数,耦合常数Jc,开始温度,结束温度,温度间隔,pi
    n=10;
    Jc=1.0;
    T_initial=0.1d0;
    T_final=2.0d0;
    T_interval=0.1d0;
    pi=3.1415926d0
    !读取:弛豫步数,采样数,间隔数
    Thermal=10000;
    sampling=1000;
    interval=1000;
    !计算温度样本数
    T_number=((T_final-T_initial)/T_interval);
    !初始温度
    T=T_initial;

    Nsite=n*n;
	
	arm=0              ! average relative magnetization
    asm=0              ! average square magnetization
    ave=0              ! average relative energy
    ase=0              ! average square energy
	!aro=0              !平均超流密度
    M_2=0              !update9.15
    M_4=0              !update9.15

    !预分配
    call input

    call neighbor

    call init_random_seed()

    do i=1,T_number
        beta=1.0/T;
        call generate_ising
        call generate_xy
        do ii=1,Thermal!弛豫
            call flip_spin
			!call flip_spin_xy  !update 9.26
        enddo
        do ii=1,sampling!采样
		sqm=rem*rem            ! calculate square magnetization
        arm=arm+rem            ! accumulate relative magnetization
        asm=asm+sqm            ! accumulate square magnetization
        ave=ave+energy            ! accumulate energy density
        ase=ase+energy*energy        ! accumulate square energy
        M_2=M_2+rem_2*rem_2    ! update 9.15
        M_4=M_4+rem_2**4       ! update 9.15
            do iii=1,interval
                call flip_spin
				!call flip_spin_xy  !update 9.26
            enddo
            call calculate
			!call calculate_xy       !update9.26
            energy_sampling(ii)=energy;
			var_sampling(ii)=var
			!rho_s_sampling(ii)=rho_s   !update9.26
        enddo
		arm=arm/sampling
        asm=asm/sampling
        ave=ave/sampling
        ase=ase/sampling
        M_2=M_2/sampling           !update9.15
        M_4=M_4/sampling           !update9.15
        susp=asm-arm*arm
        speh=(ase-ave*ave)/T/ T 
        tbr=(M_2**2)/M_4       !update9.15  
        energy=sum(energy_sampling)/sampling
		var=-1.0/2*energy-beta*sum(var_sampling)/sampling
		!write(*,*) T,energy,var
		!pause
		!rho_s=sum(rho_s_sampling)/sampling   !update9.26
        open( 1, file = 'E.dat',status="unknown",access="append")
		print'(9f12.5)',T,energy,arm,asm,susp,ase,ave,speh,tbr,var
        write(1,*) T,energy,arm,asm,susp,ase,ave,speh,tbr,var
        close(1)
        write(*,*) T,energy,arm,asm,susp,ase,ave,speh,tbr,var
        T=T+T_interval;
    enddo
    end program mcmc_mix_metropolis    